#include <algorithm>
#include <iostream>
#include <cmath>
#include <vector>
#include <stdexcept> 

using namespace std;

// 定义ppForm类，用于样条插值
class ppForm {
private:
    int numPoints;  // 插值点数量
    int boundaryConditionType;  // 边界条件类型：1=完全三次样条，2=D2三次样条，3=自然三次样条
    double boundaryCondition1, boundaryCondition2;  // 边界条件值
    vector<double> xValues, yValues;  // x和y的数据点

public:
    // 构造函数，接受x和y的数据点
    ppForm(vector<double> xData, vector<double> yData) {
        xValues = xData;
        yValues = yData;
        numPoints = xData.size();

        // 检查x和y的数据点数量是否匹配
        if (xData.size() != yData.size()) {
            throw std::invalid_argument("Error: x and y data points do not match!");
        }
    }


    // S1方法，执行线性插值
    double S1(double xx) {
        // 遍历所有数据点
        for (int i = 0; i < numPoints - 1; i++) {
            // 检查xx是否在当前的x值区间内
            if (xx >= xValues[i] && xx <= xValues[i + 1]) {
                double slope, intercept, interpolatedY;
                
                // 计算斜率和截距
                slope = (yValues[i + 1] - yValues[i]) / (xValues[i + 1] - xValues[i]);
                intercept = yValues[i] - slope * xValues[i];
                
                // 计算并返回插值结果
                interpolatedY = intercept + slope * xx;
                return interpolatedY;
            }       
            else{
                break;
            }
        
        }
    }

    
    double S3(int boundaryConditionType, double boundaryCondition1, double boundaryCondition2, double interpolationPoint) 
    {
        double coefficients[10000][4]; // 存储插值系数
        this->boundaryConditionType = boundaryConditionType;
        this->boundaryCondition1 = boundaryCondition1;
        this->boundaryCondition2 = boundaryCondition2;

        vector<double> muValues(numPoints), lambdaValues(numPoints); // 存储mu和lambda值
        // 计算mu_i和lambda_i的值
        for (int i = 1; i < numPoints - 1; i++) {
            muValues[i] = (xValues[i] - xValues[i - 1]) / (xValues[i + 1] - xValues[i - 1]);
            lambdaValues[i] = (xValues[i + 1] - xValues[i]) / (xValues[i + 1] - xValues[i - 1]);
        }

        // 差商表
        double dividedDifferenceTable[numPoints + 2][3];
        double extendedX[numPoints + 2]; // 扩展的X值数组
        dividedDifferenceTable[0][0] = yValues[0];
        dividedDifferenceTable[numPoints + 1][0] = yValues[numPoints - 1];
        extendedX[0] = xValues[0];
        extendedX[numPoints + 1] = xValues[numPoints - 1];
        for (int i = 1; i <= numPoints; i++) {
            dividedDifferenceTable[i][0] = yValues[i - 1];
            extendedX[i] = xValues[i - 1];
        }

        // 填充差商表的其他元素
        int diagonalIndex = 1; // 对角线索引
        for (int col = 1; col < 3; col++) {
            for (int row = col; row < numPoints + 2; row++) {
                if (row == 1 && col == 1) {
                    dividedDifferenceTable[row][col] = boundaryCondition1;
                } else if (row == numPoints + 1 && col == 1) {
                    dividedDifferenceTable[row][col] = boundaryCondition2;
                } else {
                    dividedDifferenceTable[row][col] = (dividedDifferenceTable[row][col - 1] - dividedDifferenceTable[row - 1][col - 1]) / (extendedX[row] - extendedX[row - diagonalIndex]);
                }
            }
            diagonalIndex += 1;
        }
        vector<double> secondDerivatives(numPoints); // 存储二阶导数值
        vector<double> rhsValues(numPoints + 1); // 右侧值

        // boundaryConditionType=1, 完全三次样条
        if (boundaryConditionType == 1) {
            for (int i = 0; i < numPoints; i++) {
                rhsValues[i + 1] = dividedDifferenceTable[i + 2][2] * 6;
            }

            // 追赶法解三对角线性方程组
            // 构造三对角矩阵
            double diagonal[numPoints + 1], lower[numPoints + 1], upper[numPoints + 1], tempY[numPoints + 1];
            for (int i = 1; i <= numPoints; i++) {
                diagonal[i] = 2;
            }
            for (int i = 1; i <= numPoints - 1; i++) {
                if (i == 1) {
                    lower[i] = 1;
                } else {
                    lower[i] = lambdaValues[i - 1];
                }
            }
            lower[numPoints] = 0;
            upper[1] = 0;
            for (int i = 2; i <= numPoints; i++) {
                if (i == numPoints) {
                    upper[i] = 1;
                } else {
                    upper[i] = muValues[i - 1];
                }
            }
            // 解线性方程组
            double alpha[numPoints + 1], beta[numPoints + 1];
            alpha[1] = diagonal[1];
            for (int i = 1; i <= numPoints - 1; i++) {
                beta[i] = lower[i] / alpha[i];
                alpha[i + 1] = diagonal[i + 1] - upper[i + 1] * beta[i];
            }
            tempY[1] = rhsValues[1] / alpha[1];
            for (int i = 2; i <= numPoints; i++) {
                tempY[i] = (rhsValues[i] - upper[i] * tempY[i - 1]) / alpha[i];
            }
            secondDerivatives[numPoints - 1] = tempY[numPoints];
            for (int i = numPoints - 1; i >= 1; i--) {
                secondDerivatives[i - 1] = tempY[i] - beta[i] * secondDerivatives[i];
            }
            coefficients[0][3] = (secondDerivatives[1] - secondDerivatives[0]) / (6 * (xValues[1] - xValues[0]));
            coefficients[0][2] = secondDerivatives[0] / 2;
            coefficients[0][1] = boundaryCondition1;
            coefficients[0][0] = yValues[0];
            for (int i = 1; i < numPoints - 1; i++) {
                coefficients[i][3] = (secondDerivatives[i + 1] - secondDerivatives[i]) / (6 * (xValues[i + 1] - xValues[i]));
                coefficients[i][2] = secondDerivatives[i] / 2;
                coefficients[i][1] = dividedDifferenceTable[i + 2][1] - (secondDerivatives[i + 1] + 2 * secondDerivatives[i]) * (xValues[i + 1] - xValues[i]) / 6;
                coefficients[i][0] = yValues[i];
            }
        }
        // boundaryConditionType=2, D2三次样条; boundaryConditionType=3, 自然三次样条
        else if (boundaryConditionType == 2 || boundaryConditionType == 3) {
            if (boundaryConditionType == 3) {
                if (boundaryCondition1 != 0 || boundaryCondition2 != 0) {
                    std::cout << "Error: Not a natural spline!" << std::endl;
                    exit(1);
                }
            }
        
        rhsValues[2] = dividedDifferenceTable[3][2] * 6 - boundaryCondition1 * muValues[1];
        rhsValues[numPoints - 1] = dividedDifferenceTable[numPoints][2] * 6 - boundaryCondition2 * lambdaValues[numPoints - 2];
        for(int i = 2; i < numPoints - 1; i++) {
            rhsValues[i + 1] = dividedDifferenceTable[i + 2][2] * 6;
        }

        double diagonal[numPoints + 1], lower[numPoints + 1], upper[numPoints + 1], solutionY[numPoints + 1];
        for(int i = 1; i <= numPoints; i++) {
            diagonal[i] = 2; // 对角线元素初始化为2
        }
        for(int i = 1; i <= numPoints - 1; i++) {
            if(i == 1) {
                lower[i] = 1; // 第一个下对角线元素为1
            } else {
                lower[i] = lambdaValues[i - 1]; // 填充下对角线元素
            }
        }
        lower[numPoints] = 0; // 最后一个下对角线元素为0
        upper[1] = 0; // 第一个上对角线元素为0
        for(int i = 2; i <= numPoints; i++) {
            if(i == numPoints) {
                upper[i] = 1; // 最后一个上对角线元素为1
            } else {
                upper[i] = muValues[i - 1]; // 填充上对角线元素
            }
        }
        // 解线性方程组
        double alpha[numPoints + 1], beta[numPoints + 1];
        alpha[1] = diagonal[1];
        for(int i = 1; i <= numPoints - 1; i++) {
            beta[i] = lower[i] / alpha[i];
            alpha[i + 1] = diagonal[i + 1] - upper[i + 1] * beta[i];
        }
        solutionY[1] = rhsValues[1] / alpha[1];
        for(int i = 2; i <= numPoints; i++) {
            solutionY[i] = (rhsValues[i] - upper[i] * solutionY[i - 1]) / alpha[i];
        }
        secondDerivatives[numPoints - 1] = solutionY[numPoints];
        for(int i = numPoints - 1; i >= 1; i--) {
            secondDerivatives[i - 1] = solutionY[i] - beta[i] * secondDerivatives[i];
        }
        coefficients[0][3] = (secondDerivatives[1] - secondDerivatives[0]) / (6 * (xValues[1] - xValues[0]));
        coefficients[0][2] = secondDerivatives[0] / 2;
        coefficients[0][1] = boundaryCondition1;
        coefficients[0][0] = yValues[0];
        for(int i = 1; i < numPoints - 1; i++) {
            coefficients[i][3] = (secondDerivatives[i + 1] - secondDerivatives[i]) / (6 * (xValues[i + 1] - xValues[i]));
            coefficients[i][2] = secondDerivatives[i] / 2;
            coefficients[i][1] = dividedDifferenceTable[i + 2][1] - (secondDerivatives[i + 1] + 2 * secondDerivatives[i]) * (xValues[i + 1] - xValues[i]) / 6;
            coefficients[i][0] = yValues[i];
        }
        }
        else if(boundaryConditionType == 2 || boundaryConditionType == 3) {
            if(boundaryConditionType == 3) {
                if(boundaryCondition1 != 0 || boundaryCondition2 != 0) {
                    std::cout << "Error: Not a natural spline!" << std::endl;
                    exit(1);
                }
            }

        rhsValues[2] = dividedDifferenceTable[3][2] * 6 - boundaryCondition1 * muValues[1];
        rhsValues[numPoints - 1] = dividedDifferenceTable[numPoints][2] * 6 - boundaryCondition2 * lambdaValues[numPoints - 2];
        for(int i = 2; i < numPoints - 1; i++) {
            rhsValues[i + 1] = dividedDifferenceTable[i + 2][2] * 6;
        }

        // 追赶法解三对角线性方程组
        // 构造三对角矩阵
        double diagonal[numPoints + 1], lower[numPoints + 1], upper[numPoints + 1], tempY[numPoints + 1];
        for(int i = 2; i <= numPoints - 1; i++) {
            diagonal[i] = 2;
        }
        for(int i = 2; i <= numPoints - 2; i++) {
            lower[i] = lambdaValues[i - 1];
        }
        for(int i = 3; i <= numPoints - 1; i++) {
            upper[i] = muValues[i - 1];
        }
        // 解线性方程组
        double alpha[numPoints + 1], beta[numPoints + 1];
        alpha[2] = diagonal[2];
        for(int i = 2; i <= numPoints - 2; i++) {
            beta[i] = lower[i] / alpha[i];
            alpha[i + 1] = diagonal[i + 1] - upper[i + 1] * beta[i];
        }
        tempY[2] = rhsValues[2] / alpha[2];
        for(int i = 3; i <= numPoints - 1; i++) {
            tempY[i] = (rhsValues[i] - upper[i] * tempY[i - 1]) / alpha[i];
        }
        secondDerivatives[0] = boundaryCondition1;
        secondDerivatives[numPoints - 2] = tempY[numPoints - 1];
        secondDerivatives[numPoints - 1] = boundaryCondition2;
        for(int i = numPoints - 2; i >= 2; i--) {
            secondDerivatives[i - 1] = tempY[i] - beta[i] * secondDerivatives[i];
        }
        for(int i = 0; i < numPoints - 1; i++) {
            coefficients[i][3] = (secondDerivatives[i + 1] - secondDerivatives[i]) / (6 * (xValues[i + 1] - xValues[i]));
            coefficients[i][2] = secondDerivatives[i] / 2;
            coefficients[i][1] = dividedDifferenceTable[i + 2][1] - (secondDerivatives[i + 1] + 2 * secondDerivatives[i]) * (xValues[i + 1] - xValues[i]) / 6;
            coefficients[i][0] = yValues[i];
        }
        }
        else {
            std::cout << "Error: Incorrect boundary condition type!" << std::endl;
            exit(1);
        }
        for(int i = 0; i <= numPoints - 1; i++) {
        if(interpolationPoint >= xValues[i] && interpolationPoint <= xValues[i + 1]) {
            // 计算插值结果
            double interpolatedY;
            interpolatedY = coefficients[i][0] + coefficients[i][1] * (interpolationPoint - xValues[i]) + 
                            coefficients[i][2] * pow(interpolationPoint - xValues[i], 2) + 
                            coefficients[i][3] * pow(interpolationPoint - xValues[i], 3);
            return interpolatedY; // 返回插值点的函数值
        }
        else{
            break;
        }
        }
    }
};

class Bspline {
private:
    int numPoints; // 插值点数量
    int boundaryConditionType; // 边界条件类型：1=完全三次样条，2=D2三次样条，3=自然三次样条
    double boundaryCondition1, boundaryCondition2; // 边界条件值
    vector<double> xValues, yValues; // x和y的数据点

public:
    // 构造函数，接受x和y的数据点
    Bspline(vector<double> _x, vector<double> _y) {
        xValues = _x;
        yValues = _y;
        numPoints = _x.size();
        if (_x.size() != _y.size()) {
            throw std::invalid_argument("Error: The sizes of x and y data points do not match!");
        }
    }

    // 计算B样条基函数
    double B(int i, int degree, double xx) {
        if (degree == 0) {
            if (xx >= xValues[i - 1] && xx <= xValues[i]) {
                return 1.0;
            } else {
                return 0;
            }
        } else {
            double term1 = (xx - xValues[i - 1]) / (xValues[i + degree - 1] - xValues[i - 1]) * B(i, degree - 1, xx);
            double term2 = (xValues[i + degree] - xx) / (xValues[i + degree] - xValues[i]) * B(i + 1, degree - 1, xx);
            return term1 + term2;
        }
    }

    // 一阶B样条插值
    double S1(double xx) {
        double result = 0.0;
        for (int i = 0; i < numPoints; i++) {
            result += yValues[i] * B(i + 1, 1, xx);
        }
        return result;
    }

    double S2(double boundaryCondition1, double boundaryCondition2, double interpolationPoint) {
    vector<double> coefficients(numPoints + 2, 0); // 初始化系数为0
    this->boundaryCondition1 = boundaryCondition1;
    this->boundaryCondition2 = boundaryCondition2;

    vector<double> b(numPoints + 1), a(numPoints + 1, 6), c(numPoints + 1, 1), d(numPoints + 1, 1), yy(numPoints + 1);
    
    // 初始化b数组
    b[1] = yValues[0] * 8 - 2 * boundaryCondition1;
    b[numPoints] = yValues[numPoints - 1] * 8 - 2 * boundaryCondition2;
    for (int i = 2; i < numPoints; i++) {
        b[i] = 8 * yValues[i - 1];
    }

    // 修改a数组的特定值
    a[1] = 5;
    a[numPoints] = 5;

    // 初始化c和d数组
    c[numPoints] = 0;
    d[1] = 0;

    // 追赶法求解线性方程组
    vector<double> alpha(numPoints + 1), beta(numPoints + 1);
    alpha[1] = a[1];
    for (int i = 1; i < numPoints; i++) {
        beta[i] = c[i] / alpha[i];
        alpha[i + 1] = a[i + 1] - d[i + 1] * beta[i];
    }
    yy[1] = b[1] / alpha[1];
    for (int i = 2; i <= numPoints; i++) {
        yy[i] = (b[i] - d[i] * yy[i - 1]) / alpha[i];
    }

    // 计算系数
    coefficients[numPoints] = yy[numPoints];
    for (int i = numPoints - 1; i >= 1; i--) {
        coefficients[i] = yy[i] - beta[i] * coefficients[i + 1];
    }
    coefficients[0] = 2 * boundaryCondition1 - coefficients[1];
    coefficients[numPoints + 1] = 2 * boundaryCondition2 - coefficients[numPoints];

    // 计算插值结果
    double result = 0;
    for (int i = 0; i < numPoints + 1; i++) {
        result += coefficients[i] * B(i, 2, interpolationPoint);
    }
    return result;
    }

    double S3(int conditionType, double boundaryCondition1, double boundaryCondition2, double interpolationPoint) {
    vector<double> coefficients(numPoints + 2, 0); // 初始化系数数组
    this->boundaryConditionType = conditionType;
    this->boundaryCondition1 = boundaryCondition1;
    this->boundaryCondition2 = boundaryCondition2;

    // 为追赶法准备数组
    vector<double> b(numPoints + 1), a(numPoints + 1), c(numPoints + 1, 1), d(numPoints + 1, 1), yy(numPoints + 1);

    if (conditionType == 1) { // 完全三次样条
        for (int i = 1; i <= numPoints - 1; i++) {
            b[i] = 6 * yValues[i - 1];
        }
        b[1] = yValues[0] * 3 + boundaryCondition1;
        b[numPoints] = yValues[numPoints - 1] * 3 - boundaryCondition2;
        for (int i = 1; i <= numPoints - 1; i++) {
            a[i] = 4;
        }
        a[1] = 2;
        a[numPoints] = 2;
    } else if (conditionType == 2 || conditionType == 3) { // D2三次样条或自然三次样条
        if (conditionType == 3 && (boundaryCondition1 != 0 || boundaryCondition2 != 0)) {
            throw std::invalid_argument("Error: Not a natural spline!");
        }
        b[2] = 6 * yValues[1] - boundaryCondition1;
        b[numPoints - 1] = 6 * yValues[numPoints - 2] - boundaryCondition2;
        for (int i = 3; i <= numPoints - 2; i++) {
            b[i] = 6 * yValues[i - 1];
        }
        for (int i = 2; i <= numPoints - 1; i++) {
            a[i] = 4;
        }
        a[1] = 1;
        a[numPoints] = 1;
    } else {
        throw std::invalid_argument("Error: Incorrect boundary condition type!");
    }

    // 追赶法求解线性方程组
    vector<double> alpha(numPoints + 1), beta(numPoints + 1);
    alpha[1] = a[1];
    for (int i = 1; i < numPoints; i++) {
        beta[i] = c[i] / alpha[i];
        alpha[i + 1] = a[i + 1] - d[i + 1] * beta[i];
    }
    yy[1] = b[1] / alpha[1];
    for (int i = 2; i <= numPoints; i++) {
        yy[i] = (b[i] - d[i] * yy[i - 1]) / alpha[i];
    }
    coefficients[numPoints] = yy[numPoints];
    for (int i = numPoints - 1; i >= 1; i--) {
        coefficients[i] = yy[i] - beta[i] * coefficients[i + 1];
    }
    coefficients[0] = 2 * boundaryCondition1 - coefficients[1];
    coefficients[numPoints + 1] = 2 * boundaryCondition2 - coefficients[numPoints];

    // 计算插值结果
    double result = 0;
    for (int i = 0; i < numPoints + 2; i++) {
        result += coefficients[i] * B(i - 1, 3, interpolationPoint);
    }
    return result;
}
};

    